
#### HOUSEKEEPING ####

rm(list = ls())

package_list <- c("data.table", "tidyverse", "latex2exp")
lapply(package_list, require, character.only = TRUE)

nsim       <- 200
just_graph <- 1
bw         <- 2

source("00_paths.R")

## MASKED PARAMETERS
n_tier <- NA  # Number of service tiers
min_tier <- NA  # Minimum tier number
max_tier <- NA  # Maximum tier number

#### SET UP ####

set.seed(6392)

adaptive_kernel <- function(x,y,bw,gridx,adapt=TRUE) {
  n <- length(gridx)
  kernel <- function(u) {
    return((1/sqrt(2*pi))*exp(-(u^2)/2))
  }
  adapt_bw <- function(ig) {
    dist <- abs(x-gridx[ig])
    local_dens <- 1/mean(dist[dist<bw])
    return(bw/local_dens)
  }
  smooth <- function(ig) {
    cur_bw <- ifelse(adapt, adapt_bw(ig), bw)
    w <- kernel((x-gridx[ig])/cur_bw)
    return(sum(w*y)/sum(w))
  }
  ys <- sapply(1:n, smooth)
}

## LOOP OVER SIMULATIONS

if (just_graph==0) {
    
    for (task_id in c(0:nsim)) {
      
      print(task_id)
      
      ## LOAD CURRENT SIM RESULTS
      file_list <- file.path(res_dir, paste0(sprintf("%03d", task_id), "-cf.csv"))
      df        <- do.call(rbind, lapply(file_list, fread))
      
      ## RESAMPLE TREATED HHS
      if (task_id>0) {
        nt    <- nrow(df)
        t_samp <- sample(1:nt, nt, replace = TRUE)
        df <- df[t_samp,]
      }
      
      ## GRAPHS
      x_var    <- "e_ovr"
      df$x_var <- df[[x_var]]
      
      tgraph1_df <- df %>%
        mutate(bill_i = (en_bill_i - st_bill_i - (cf_en_bill_i - cf_st_bill_i))) %>%
        mutate(bill_i0 = ((cf_en_bill_i - cf_st_bill_i))) %>%
        mutate(bill_t = (en_bill_t - st_bill_t - (cf_en_bill_t - cf_st_bill_t))) %>%
        mutate(bill_t0 = ((cf_en_bill_t - cf_st_bill_t))) %>%
        mutate(bill_u = (en_ovr)) %>%
        mutate(bill_u0 = (0)) %>%
        mutate(bill = bill_i + bill_t + bill_u) %>%
        mutate(bill0 = bill_i0 + bill_t0 + bill_u0) %>%
        mutate(tieru = upg_tier - cf_upg_tier) %>%
        mutate(tieru0 = cf_upg_tier) %>%
        mutate(tieru = ifelse(st_tier<max_tier, tieru, NA)) %>%
        mutate(tieru0 = ifelse(st_tier<max_tier, tieru0, NA)) %>%
        mutate(tierd = dng_tier - cf_dng_tier) %>%
        mutate(tierd0 = cf_dng_tier) %>%
        mutate(tierd = ifelse(st_tier>min_tier, tierd, NA)) %>%
        mutate(tierd0 = ifelse(st_tier>min_tier, tierd0, NA)) %>%
        mutate(tvd = drop_vid - cf_drop_vid) %>%
        mutate(tvd0 = cf_drop_vid) %>%
        mutate(tvd = ifelse(st_vid==1, tvd, NA)) %>%
        mutate(tvd0 = ifelse(st_vid==1, tvd0, NA)) %>%
        mutate(tva = add_vid - cf_add_vid) %>%
        mutate(tva0 = cf_add_vid) %>%
        mutate(tva = ifelse(st_vid==0, tva, NA)) %>%
        mutate(tva0 = ifelse(st_vid==0, tva0, NA)) %>%
        mutate(dgb = (en_tot_gb - cf_en_gb)) %>%
        mutate(dgb0 = (cf_en_gb)) %>%
        mutate(dgb_video = (en_gb_video - cf_en_video)) %>%
        mutate(dgb_browsing = (en_gb_browsing - cf_en_browsing)) %>%
        mutate(dgb_other = (en_gb_other - cf_en_other)) %>%
        mutate(dgb_netflix = (en_gb_netflix - cf_en_netflix)) %>%
        mutate(dgb_youtube = (en_gb_youtube - cf_en_youtube)) %>%
        mutate(dgb_linear = (en_gb_hulu + en_gb_slingtv - cf_en_hulu - cf_en_slingtv)) %>%
        mutate(dgb_othvid = (dgb_video - dgb_netflix - dgb_youtube - dgb_linear)) %>%
        dplyr::select(x=e_ovr, bill, bill0, bill_i, bill_i0, bill_t, bill_t0, bill_u, bill_u0,
                      starts_with("tieru"), starts_with("tierd"),
                      starts_with("tvd"), starts_with("tva"), starts_with("dgb"))
      
        var_list <- setdiff(names(tgraph1_df), "x")
        nv <- length(var_list)
        ng <- 26
        gridx <- seq(0,50,length.out=ng)
        sim1_df <- data.frame(
          v = rep(var_list, each=ng),
          s = rep(task_id, times=ng*nv),
          x = rep(gridx, times=nv),
          y = NA
        )
        
        for (cur_v in var_list) {
          
          x <- tgraph1_df$x
          y <- tgraph1_df[[cur_v]]
          
          keep <- !is.na(y)
          x <- x[keep]
          y <- y[keep]
          
          sim1_df$y[sim1_df$v==cur_v] <- adaptive_kernel(x,y,bw=bw,gridx,adapt=T)
          
        }
      
      
      tgraph2_df <- df %>%
        mutate(dgb = (en_tot_gb - cf_en_gb)) %>%
        mutate(dgb0 = (cf_en_gb)) %>%
        mutate(dgb_video = (en_gb_video - cf_en_video)) %>%
        mutate(dgb_browsing = (en_gb_browsing - cf_en_browsing)) %>%
        mutate(dgb_other = (en_gb_other - cf_en_other)) %>%
        mutate(dgb_netflix = (en_gb_netflix - cf_en_netflix)) %>%
        mutate(dgb_youtube = (en_gb_youtube - cf_en_youtube)) %>%
        mutate(dgb_linear = (en_gb_hulu + en_gb_slingtv - cf_en_hulu - cf_en_slingtv)) %>%
        mutate(dgb_othvid = (dgb_video - dgb_netflix - dgb_youtube - dgb_linear)) %>%
        mutate(dgb_video0 = (cf_en_video)) %>%
        mutate(dgb_browsing0 = (cf_en_browsing)) %>%
        mutate(dgb_other0 = (cf_en_other)) %>%
        mutate(dgb_netflix0 = (cf_en_netflix)) %>%
        mutate(dgb_youtube0 = (cf_en_youtube)) %>%
        mutate(dgb_linear0 = (cf_en_hulu + cf_en_slingtv)) %>%
        mutate(dgb_othvid0 = (dgb_video0 - dgb_netflix0 - dgb_youtube0 - dgb_linear0)) %>%
        dplyr::select(x=x_var, upg=upg_tier, starts_with("dgb"))
      
      var_list <- setdiff(names(tgraph2_df), c("x", "upg"))
      nv <- length(var_list)
      ng <- 26
      gridx <- seq(0,50,length.out=ng)
      tsim_df <- data.frame(
        v = rep(var_list, each=ng),
        s = rep(task_id, times=ng*nv),
        x = rep(gridx, times=nv),
        y = NA
      )
      sim2_df <- rbind(tsim_df %>% mutate(u=1), tsim_df %>% mutate(u=0))
      
      for (cur_v in var_list) {
        
        for (cur_u in c(0,1)) {
          
          x <- tgraph2_df$x[tgraph2_df$upg==cur_u]
          y <- tgraph2_df[[cur_v]][tgraph2_df$upg==cur_u]
          keep <- !is.na(y)
          x <- x[keep]
          y <- y[keep]
          
          sim2_df$y[sim2_df$v==cur_v & sim2_df$u==cur_u] <- adaptive_kernel(x,y,bw=bw,gridx,adapt=T)
          
        }
        
      }
      
      if (!exists("graph1_df")) {
        graph1_df <- sim1_df
      } else {
        graph1_df <- graph1_df %>% rbind(sim1_df)
      }
      
      if (!exists("graph2_df")) {
        graph2_df <- sim2_df
      } else {
        graph2_df <- graph2_df %>% rbind(sim2_df)
      }
      

      ### COMPUTE ATES FOR EACH RESAMPLE

      # mean effects
      ate_df <- df %>%
        mutate(bill_i = (en_bill_i - st_bill_i - (cf_en_bill_i - cf_st_bill_i))) %>%
        mutate(bill_i0 = ((cf_en_bill_i - cf_st_bill_i))) %>%
        mutate(bill_t = (en_bill_t - st_bill_t - (cf_en_bill_t - cf_st_bill_t))) %>%
        mutate(bill_t0 = ((cf_en_bill_t - cf_st_bill_t))) %>%
        mutate(bill_u = (en_ovr)) %>%
        mutate(bill_u0 = (0)) %>%
        mutate(bill = bill_i + bill_t + bill_u) %>%
        mutate(bill0 = bill_i0 + bill_t0 + bill_u0) %>%
        mutate(tieru = upg_tier - cf_upg_tier) %>%
        mutate(tieru0 = cf_upg_tier) %>%
        mutate(tieru = ifelse(st_tier<max_tier, tieru, NA)) %>%
        mutate(tieru0 = ifelse(st_tier<max_tier, tieru0, NA)) %>%
        mutate(tierd = dng_tier - cf_dng_tier) %>%
        mutate(tierd0 = cf_dng_tier) %>%
        mutate(tierd = ifelse(st_tier>min_tier, tierd, NA)) %>%
        mutate(tierd0 = ifelse(st_tier>min_tier, tierd0, NA)) %>%
        mutate(tvd = drop_vid - cf_drop_vid) %>%
        mutate(tvd0 = cf_drop_vid) %>%
        mutate(tvd = ifelse(st_vid==1, tvd, NA)) %>%
        mutate(tvd0 = ifelse(st_vid==1, tvd0, NA)) %>%
        mutate(tva = add_vid - cf_add_vid) %>%
        mutate(tva0 = cf_add_vid) %>%
        mutate(tva = ifelse(st_vid==0, tva, NA)) %>%
        mutate(tva0 = ifelse(st_vid==0, tva0, NA)) %>%
        mutate(dgb = (en_tot_gb - cf_en_gb)) %>%
        mutate(dgb0 = (cf_en_gb)) %>%
        mutate(dgb_video = (en_gb_video - cf_en_video)) %>%
        mutate(dgb_browsing = (en_gb_browsing - cf_en_browsing)) %>%
        mutate(dgb_other = (en_gb_other - cf_en_other)) %>%
        mutate(dgb_netflix = (en_gb_netflix - cf_en_netflix)) %>%
        mutate(dgb_youtube = (en_gb_youtube - cf_en_youtube)) %>%
        mutate(dgb_linear = (en_gb_hulu + en_gb_slingtv - cf_en_hulu - cf_en_slingtv)) %>%
        mutate(dgb_othvid = (dgb_video - dgb_netflix - dgb_youtube - dgb_linear)) %>%
        mutate(dgb_video0 = (cf_en_video)) %>%
        mutate(dgb_browsing0 = (cf_en_browsing)) %>%
        mutate(dgb_other0 = (cf_en_other)) %>%
        mutate(dgb_netflix0 = (cf_en_netflix)) %>%
        mutate(dgb_youtube0 = (cf_en_youtube)) %>%
        mutate(dgb_linear0 = (cf_en_hulu + cf_en_slingtv)) %>%
        mutate(dgb_othvid0 = (dgb_video0 - dgb_netflix0 - dgb_youtube0 - dgb_linear0))

      cate_df <- ate_df

      ate_df <- ate_df %>%
        dplyr::select(x=e_ovr, bill, bill0, starts_with("bill_i"), starts_with("bill_t"),
                      starts_with("bill_u"), starts_with("tieru"), starts_with("tierd"),
                      starts_with("tvd"), starts_with("tva"), starts_with("dgb")) %>%
        select(-x) %>%
        summarise(across(everything(), mean, na.rm=TRUE)) %>%
        t() %>%
        as.data.frame()

      ate_df$Var <- rownames(ate_df)
      ate_df$ATE <- ifelse(grepl("0", ate_df$Var), NA, ate_df$V1)
      ate_df$CF <- ifelse(grepl("0", ate_df$Var), ate_df$V1, NA)
      ate_df$Var <- gsub("0", "", ate_df$Var)

      ate_df <- ate_df %>%
        group_by(Var) %>%
        mutate(ATE = max(ATE, na.rm=TRUE)) %>%
        mutate(CF = max(CF, na.rm=TRUE)) %>%
        ungroup() %>%
        select(Var, ATE, CF) %>%
        distinct() %>%
        mutate(Sim = task_id)

      cate_df <- cate_df %>%
        dplyr::select(x=e_ovr, bill, bill0, starts_with("bill_i"), starts_with("bill_t"),
                      starts_with("bill_u"), starts_with("tieru"), starts_with("tierd"),
                      starts_with("tvd"), starts_with("tva"), starts_with("dgb"), st_vid, en_vid) %>%
        select(-x) %>%
        group_by(st_vid, en_vid) %>%
        mutate(N = n()) %>%
        summarise(across(everything(), mean, na.rm=TRUE)) %>%
        t() %>%
        as.data.frame()

      cate_df$Var <- rownames(cate_df)
      cate_df$ATE <- ifelse(grepl("0", cate_df$Var), NA, cate_df$V1)
      cate_df$CF <- ifelse(grepl("0", cate_df$Var), cate_df$V1, NA)
      cate_df$Var <- gsub("0", "", cate_df$Var)


      if (!exists("graph4_df")) {
        graph4_df <- ate_df
      } else {
        graph4_df <- graph4_df %>% rbind(ate_df)
      }

      fwrite(x = graph4_df, file = file.path(dat_dir, "ate_resamp.csv"))

    fwrite(x = graph1_df, file = file.path(dat_dir, paste0("smoothgraph1.csv")))
    fwrite(x = graph2_df, file = file.path(dat_dir, paste0("smoothgraph2.csv")))

  rm(graph1_df)
  rm(graph2_df)
  
  }

} else {

  graph1_df <- fread(file.path(dat_dir, paste0("smoothgraph1_bw300.csv")))
  graph2_df <- fread(file.path(dat_dir, paste0("smoothgraph2_bw300.csv")))
  ate_df <- fread(file.path(dat_dir, "ate_resamp.csv"))

}
  

#### FIGURES ####

## GRAPHING PARAMS
cpalette3 <- c('gray20','gray50','gray80')
lt3       <- c("solid","dashed","dotted")
markers3  <- c(16,15,17)
cpalette4 <- c('gray20','gray40','gray60','gray80')
lt4       <- c("solid","dashed","dotted","dotdash")
markers4  <- c(16,15,17,18)
  
## BASE GRAPH
base_graph <- list(geom_hline(yintercept = 0, color="gray90"),
                   geom_line(aes(y = y,
                                  x = x)),
                   geom_ribbon(aes(x = x,
                                     ymin = lb,
                                     ymax = ub),
                                 fill="gray50", alpha=0.2),
                   theme_bw(),
                   theme(axis.text=element_text(size=rel(0.75)),
                         axis.title=element_text(size=rel(0.75)),
                         legend.position = c(0.16,0.86),
                         legend.text=element_text(size=rel(0.6)),
                         text = element_text(family = "Times New Roman"),
                         panel.grid.minor = element_blank()
                   ))

#### ~ FIGURES 8a, 4a, 4b, 4c, 4d, 5 ####

dv_df <- data.frame(
  dv   = c("tieru", "tierd", "tvd", "tva", "dgb", "bill"),
  ylab = c("Tier Upgrade Probability", "Tier Downgrade Probability", "TV Drop Probability", 
           "TV Add Probability", "Usage (GB)", "Monthly Bill ($)"),
  fname = c("Figure4c", "Figure4d", "Figure4b", "Figure4a", 
            "Figure5", "Figure8a")
)

for (idv in 1:nrow(dv_df)) {
  
  graph_df <- graph1_df %>%
    filter(v==dv_df$dv[idv]) %>%
    group_by(x) %>%
    mutate(lb = quantile(y, probs=0.025)) %>%
    mutate(ub = quantile(y, probs=0.975)) %>%
    ungroup() %>%
    filter(s==0) %>%
    select(x,y,lb,ub)
  
  g1 <- ggplot(graph_df) +
    base_graph +
    labs(
      x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
      y = dv_df$ylab[idv]
    )
  
  if (dv_df$dv[idv] %in% c("tvd","tva")) {
    g1 <- g1 + coord_cartesian(ylim=c(-0.04,0.1))
  }
  
  ggsave(plot=g1,
         filename=paste0(fig_dir, "/", dv_df$fname[idv], ".png"),
         height=3, width=3, device="png")
  
}


#### ~ FIGURE 8b ####

graph_df <- graph1_df %>%
  filter(v %in% c("bill_i", "bill_t", "bill_u")) %>%
  group_by(x,v) %>%
  mutate(lb = quantile(y, probs=0.025)) %>%
  mutate(ub = quantile(y, probs=0.975)) %>%
  ungroup() %>%
  filter(s==0) %>%
  select(x,y,lb,ub,v) %>%
  mutate(type = ifelse(v=="bill_u", "Usage Fees", ifelse(v=="bill_i", "Internet Plan", "TV Plan"))) %>%
  mutate(type = factor(type, levels=c("Usage Fees","Internet Plan","TV Plan")))

ggplot(graph_df) +
  aes(x=x, y = y, color = as.factor(type), shape = as.factor(type), linetype = as.factor(type)) +
  geom_hline(yintercept = 0, color = "gray") +
  geom_line(aes(x=x,y=y)) +
  geom_ribbon(aes(ymin = lb, ymax = ub, fill=as.factor(type)), alpha=0.2, colour=NA) +
  theme_bw() +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Monthly Bill ($)",
    color = ""
  ) +
  scale_color_manual(values = cpalette3, name = "Legend", labels = levels(graph_df$type)) +
  scale_fill_manual(values = cpalette3, name = "Legend", labels = levels(graph_df$type)) +
  scale_shape_manual(values = c(16,15,17), name = "Legend", labels = levels(graph_df$type)) +
  scale_linetype_manual(values = lt3, name = "Legend", labels = levels(graph_df$type)) +
  theme(legend.position = c(0.16,0.86),
        legend.text=element_text(size=rel(0.6)),
        axis.text=element_text(size=rel(0.75)),
        axis.title=element_text(size=rel(0.75)),
        legend.title=element_blank(),
        legend.margin = margin(0, 0, 0, 0),
        legend.spacing.x = unit(0, "mm"),
        legend.spacing.y = unit(0, "mm"),
        legend.key.height = unit(3, "mm"),
        legend.key.width = unit(6, "mm"),
        text = element_text(family = "Times New Roman"),
        panel.grid.minor = element_blank())

ggsave(filename=paste0(fig_dir, "/", "Figure8b", ".png"),
       height=3, width=3, device="png")

#### ~ FIGURE 6 ####

graph_df <- graph2_df %>%
  filter(v=="dgb") %>%
  group_by(x,u) %>%
  mutate(lb = quantile(y, probs=0.025)) %>%
  mutate(ub = quantile(y, probs=0.975)) %>%
  ungroup() %>%
  filter(s==0) %>%
  select(x,y,lb,ub,u)

common_y_scale <- c(-8.25,4)
common_y_ticks <- seq(-8,4, by=2)

## FIGURE 6a

ggplot(graph_df %>% filter(u==0)) +
  base_graph +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Usage (GB)"
  ) +
  coord_cartesian(ylim = common_y_scale) +
  scale_y_continuous(breaks = common_y_ticks)

ggsave(filename=paste0(fig_dir, "/", "Figure6a", ".png"),
       height=3, width=3, device="png")

## FIGURE 6b

ggplot(graph_df %>% filter(u==1)) +
  base_graph +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Usage (GB)"
  ) +
  coord_cartesian(ylim = common_y_scale) +
  scale_y_continuous(breaks = common_y_ticks)

ggsave(filename=paste0(fig_dir, "/", "Figure6b", ".png"),
       height=3, width=3, device="png")

#### ~ FIGURE 7 ####

x_lab <- "EXPECTED OVERAGE PERCENTILE BIN"

graph_df <- graph2_df %>%
  filter(grepl("dgb_", v)) %>%
  filter(!grepl("0", v)) %>%
  group_by(x,v,u) %>%
  mutate(lb = quantile(y, probs=0.025)) %>%
  mutate(ub = quantile(y, probs=0.975)) %>%
  ungroup() %>%
  filter(s==0) %>%
  select(x,y,lb,ub,v,u) %>%
  mutate(type = ifelse(v=="dgb_linear", "Hulu/Sling TV", 
                       ifelse(v=="dgb_netflix", "Netflix", 
                              ifelse(v=="dgb_othvid", "Other Video",
                                     ifelse(v=="dgb_youtube", "YouTube",
                                            ifelse(v=="dgb_browsing", "Browsing",
                                                   ifelse(v=="dgb_video", "Video",
                                                          ifelse(v=="dgb_other", "Other", v)))))))) %>%
  mutate(type = factor(type, levels=c("Video", "Browsing", "Other", 
                                      "Netflix", "YouTube", "Hulu/Sling TV", "Other Video"))) %>%
  mutate(sig = factor(ifelse(!(0>lb & 0<ub), "Sig", "Not Sig"),
                      levels = c("Sig", "Not Sig")))
  
common_y_scale <- c(-4.5,3)
common_y_ticks <- seq(-4,4, by=2)

## FIGURE 7b

ggplot(graph_df %>% filter(type %in% c("Video","Browsing","Other")) %>% filter(u==1)) +
  aes(x=x, y = y, color = as.factor(type), shape = as.factor(type), linetype = as.factor(type)) +
  geom_hline(yintercept = 0, color = "gray") +
  geom_line(aes(x=x,y=y)) +
  geom_ribbon(aes(ymin = lb, ymax = ub, fill=as.factor(type)), alpha=0.2, colour=NA) +
  theme_bw() +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Usage (GB)",
    color = ""
  ) +
  scale_color_manual(values = cpalette3, name = "Legend", labels = levels(graph_df$type)) +
  scale_fill_manual(values = cpalette3, name = "Legend", labels = levels(graph_df$type)) +
  scale_shape_manual(values = markers3, name = "Legend", labels = levels(graph_df$type)) +
  scale_linetype_manual(values = lt3, name = "Legend", labels = levels(graph_df$type)) +
  theme(legend.position = c(0.16,0.1),
        legend.text=element_text(size=rel(0.6)),
        axis.text=element_text(size=rel(0.75)),
        axis.title=element_text(size=rel(0.75)),
        legend.title=element_blank(),
        legend.margin = margin(0, 0, 0, 0),
        legend.spacing.x = unit(0, "mm"),
        legend.spacing.y = unit(0, "mm"),
        legend.key.height = unit(3, "mm"),
        legend.key.width = unit(6, "mm"),
        text = element_text(family = "Times New Roman"),
        panel.grid.minor = element_blank())

ggsave(filename=paste0(fig_dir, "/", "Figure7b", ".png"),
       height=3, width=3, device="png")

## FIGURE 7a

ggplot(graph_df %>% filter(type %in% c("Video","Browsing","Other")) %>% filter(u==0)) +
  aes(x=x, y = y, color = as.factor(type), shape = as.factor(type), linetype = as.factor(type)) +
  geom_hline(yintercept = 0, color = "gray") +
  geom_line(aes(x=x,y=y)) +
  geom_ribbon(aes(ymin = lb, ymax = ub, fill=as.factor(type)), alpha=0.2, colour=NA) +
  theme_bw() +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Usage (GB)",
    color = ""
  ) +
  scale_color_manual(values = cpalette3, name = "Legend", labels = levels(graph_df$type)) +
  scale_fill_manual(values = cpalette3, name = "Legend", labels = levels(graph_df$type)) +
  scale_shape_manual(values = markers3, name = "Legend", labels = levels(graph_df$type)) +
  scale_linetype_manual(values = lt3, name = "Legend", labels = levels(graph_df$type)) +
  theme(legend.position = c(0.16,0.1),
        legend.text=element_text(size=rel(0.6)),
        axis.text=element_text(size=rel(0.75)),
        axis.title=element_text(size=rel(0.75)),
        legend.title=element_blank(),
        legend.margin = margin(0, 0, 0, 0),
        legend.spacing.x = unit(0, "mm"),
        legend.spacing.y = unit(0, "mm"),
        legend.key.height = unit(3, "mm"),
        legend.key.width = unit(6, "mm"),
        text = element_text(family = "Times New Roman"),
        panel.grid.minor = element_blank())

ggsave(filename=paste0(fig_dir, "/", "Figure7a", ".png"),
       height=3, width=3, device="png")

common_y_scale <- c(-2.25,3.25)
common_y_ticks <- seq(-2,3, by=1)

## FIGURE 7c

ggplot(graph_df %>% filter(!type %in% c("Video","Browsing","Other")) %>% filter(u==1)) +
  aes(x=x, y = y, color = as.factor(type), shape = as.factor(type), linetype = as.factor(type)) +
  geom_hline(yintercept = 0, color = "gray") +
  geom_line(aes(x=x,y=y)) +
  theme_bw() +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Usage (GB)",
    color = ""
  ) +
  scale_color_manual(values = cpalette4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  scale_fill_manual(values = cpalette4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  scale_shape_manual(values = markers4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  scale_linetype_manual(values = lt4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  theme(legend.position = c(0.2,0.88),
        legend.text=element_text(size=rel(0.6)),
        axis.text=element_text(size=rel(0.75)),
        axis.title=element_text(size=rel(0.75)),
        legend.title=element_blank(),
        legend.margin = margin(0, 0, 0, 0),
        legend.spacing.x = unit(0, "mm"),
        legend.spacing.y = unit(0, "mm"),
        legend.key.height = unit(3, "mm"),
        legend.key.width = unit(6, "mm"),
        text = element_text(family = "Times New Roman"),
        panel.grid.minor = element_blank())

ggsave(filename=paste0(fig_dir, "/", "Figure7c", ".png"),
       height=3, width=3, device="png")

## FIGURE 7d

ggplot(graph_df %>% filter(!type %in% c("Video","Browsing","Other")) %>% filter(u==0)) +
  aes(x=x, y = y, color = as.factor(type), shape = as.factor(type), linetype = as.factor(type)) +
  geom_hline(yintercept = 0, color = "gray") +
  geom_line(aes(x=x,y=y)) +
  theme_bw() +
  labs(
    x = TeX("Expected Overage $\\hat{\\rho}_i(\\lambda)$"),
    y = "Usage (GB)",
    color = ""
  ) +
  scale_color_manual(values = cpalette4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  scale_fill_manual(values = cpalette4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  scale_shape_manual(values = markers4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  scale_linetype_manual(values = lt4, name = "Legend", labels = levels(graph_df$type)[4:7]) +
  theme(legend.position = c(0.2,0.12),
        legend.text=element_text(size=rel(0.6)),
        axis.text=element_text(size=rel(0.75)),
        axis.title=element_text(size=rel(0.75)),
        legend.title=element_blank(),
        legend.margin = margin(0, 0, 0, 0),
        legend.spacing.x = unit(0, "mm"),
        legend.spacing.y = unit(0, "mm"),
        legend.key.height = unit(3, "mm"),
        legend.key.width = unit(6, "mm"),
        text = element_text(family = "Times New Roman"),
        panel.grid.minor = element_blank())

ggsave(filename=paste0(fig_dir, "/", "Figure7d", ".png"),
       height=3, width=3, device="png")

#### ~ APPENDIX ####

## FIGURE B.1

file_list <- file.path(res_dir, paste0(sprintf("%03d", 0), "-cf2.csv"))
df        <- do.call(rbind, lapply(file_list, fread))

df$e_ovr <- ifelse(df$e_ovr>100, 100, df$e_ovr)

ggplot(df %>% filter(e_ovr>0)) +
  aes(x = e_ovr) +
  geom_histogram(binwidth = 5) +
  theme_bw() +
  labs(
    x = TeX("$\\hat{\\rho}_i(\\lambda)$"),
    y = "Count",
  ) +
  theme(axis.text=element_text(size=rel(0.75)),
        axis.title=element_text(size=rel(0.75)),
        text = element_text(family = "Times New Roman"))

ggsave(filename=paste0(fig_dir, "/", "FigureB1", ".png"),
       height=3, width=3, device="png")
